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We present a new technique in order to quantify the dynamics of spatially extended systems. Using 
a test on the existence of unstable periodic orbits, we identify intermediate spatial scales, wherein 
the dynamics is characterized by maximum nontrivial determinism. This method is applied to 
earthquake catalogues containing time, coordinates and magnitude. As a result we extract a set of 
areas with significant deterministic and low-dimensional dynamics from the data. Finally, a simple 
model is used to show that these scales can be interpreted as local spatial coupling strengths. 
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I. INTRODUCTION 

In recent years a lot of techniques to analyse single 
complex time series have been developed (cf. (I)). But a 
special challenge is the analysis of spatiotemporal dy¬ 
namics, especially of natural systems. Although such 
systems become more and more important, e.g. in en¬ 
vironmental research J2| or brain imaging techniques || , 
little is known about how to analyse the corresponding 
data. Mostly the spatial extension is not taken into ac¬ 
count and the more-dimensional data are simplified to 
one-dimensional time series |Q. This may be senseful 
for systems, which behave approximately homogeneously 
in space, but in general this assumption is not fulfilled 
for natural systems. As a consequence much informa¬ 
tion is lost after the modification of the data, e.g. by 
describing a complex spatial pattern by a single number 
like a mean value or a variance or a more complicated 
parameter like a fractal dimension. More dynamical ap¬ 
proaches are based on a decomposition of spatiotemporal 
patterns into special basis functions, e.g. wavelet trans¬ 
formation H or the Karhunen-Loeve method jo). But 
these techniques are not appropriate in the case of rather 
irregular and noisy data. 

The purpose of this contribution is to search for charac¬ 
teristic spatial scales, on which the interesting dynamical 
properties of the system can be observed. This idea al¬ 
lows to take the spatial extension of a system into account 
and is applicable to a large variety of data. Analysing the 
dynamics on these scales will give a maximum of nonlin¬ 
ear and low-dimensional determinism jzj . As a result the 
underlying dynamics can be represented locally by a vec¬ 
tor in a low-dimensional space. This idea was suggested 
by Rand and Wilson || for systems, where the global 
dynamics is in a steady state and, therefore, trivial in 
the thermodynamic limit. Most natural systems, how¬ 
ever, are far from thermodynamic equilibrium |j] and 
their system size is far from infinite size. As a conse¬ 
quence, we observe complex dynamical behavior even on 
the largest observable scales. Therefore, the concept of 


Rand and Wilson has to be modified for the analysis of 
this class of systems. 

On small scales the dynamics is dominated by in¬ 
trinsic stochasticity and on large scales spatial averag¬ 
ing over dynamically desynchronized parts of the sys¬ 
tem suppresses determinism as well as stochastic fluc¬ 
tuations. Therefore, the averaging procedure itself may 
be exploited for the analysis of spatiotemporal time se¬ 
ries. We search for intermediate scales, where, on the one 
hand a deterministic signal is observed and on the other 
hand the loss of dynamical information arising from spa¬ 
tial averaging is as small as possible. To identify such 
a nontrivial determinism, we use a formalism of Pei and 
Moss |t]] based on the existence of unstable periodic or¬ 
bits. 

As a dynamically rich and interesting subject, we ap¬ 
ply this approach to natural seismicity which provides a 
lot of chaotic and fractal features, because the underly¬ 
ing processes like stress accumulation and ruptures are 
strongly nonlinear (cf. |lO|jl l| ] ). Here we mention the 
famous Gutenberg-Richter law jl3j 

logN = am. + b , (1) 

where N is the number of earthquakes with magnitude 
greater than or equal to m. The magnitude m is related 
to the seismic energy E by 

logE = cm + d (2) 

with constants a, b , c and d. Inserting Eq. (|j) into Eq. (i) 
a power-law for the density function D(E) is obtained 

D(E) = ^ = - 10 b ~^ E '- 1 =C E T . (3) 
dE c 

In general power-laws with noninteger exponent r indi¬ 
cate fractal distributions Jll] . 

Due to technical difficulties, the analysis of earthquake 
data is not straightforward. These data are point-like 
and not equidistant in time and space, of course. Fur¬ 
thermore most catalogues are not complete in the sense 
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that they contain all micro-quakes. In the next section 
we describe a pre-processing of the catalogue in order to 
generate regular data. In Sec. Ill we give the algorithm 
for the extraction of the scales. This algorithm is applied 
in Sec. IV to the real data and in Sec. |v| to the model 
data. Finally, we summarize the main ideas and results 
of our work. 


II. DATA AND PRE-PROCESSING 

The data, which are investigated here incorporate 
10.779 earthquakes recorded in Armenia between 1974 
and 1994 ©■ Each earthquake is described by a four¬ 
dimensional vector consisting of time, longitude, latitude 
and magnitude. The events are distributed very com¬ 
plex, i.e. not equidistant, in time and space. For the 
data analysis it is helpful to handle with data, which are 
at least equidistant in time. Therefore, we subdivide the 
time into intervals T) = i ■ At and the space into cells 
Aj where the local dynamics is considered; the shape of 
the spatial cells will be determined in the next section. 
Corresponding to Eq. (|2) we introduce the total energy 
Eij in this cell by 

E i:j = 10 m(5?,t) . (4) 

t£Ti,x£Aj 


For simplicity we have defined c = 1 and d = 0 for the 
constants in Eq. ©• Now we can define an effective mag¬ 
nitude M, j as a continuous function of time and space by 

Mij = log E^. (5) 

The subdivision into cells has to be chosen fine enough 
so that the loss of dynamical information is as small as 
possible. 

For a fixed spatial area A around a point xj = (lati¬ 
tude,longitude) we compute the function 

t 

VA,x f {t) = J - M A ^ f ) dt ( 6 ) 

o 

where Ma, s f is a sliding temporal mean of (t). We 

have found that 500 time intervals (corresponding to an 
interval length of 15 days) and a window length of 100 for 
the sliding mean are appropriate values for our calcula¬ 
tions. In Fig. [I] we show the spatial distribution of seismic 
activity E ? - integrated over the whole time [0, t ma x] 

E j = log J 10 m(s>t) dt. (7) 
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FIG. 1. Spatial distribution of seismic activity (colour 
coded). 


III. CHARACTERISTIC SPATIAL SCALES 

The challenge here is to include the spatial extension 
of the system in the analysis of the earthquake data. For 
this aim we want to extract spatial regions from the data, 
wherein the dynamics is highly deterministic. In contrast 
to global descriptions, where the spatial extension is sim¬ 
plified to a number, much less dynamical information is 
lost then. 

The main idea is to search for deterministic dynamics 
in time series <7Aj,x(t) corresponding to a fixed point x 
and different areas Aj. 

We call a scale characteristic , if the dynamics in one 
area is deterministic with higher significance than in the 
other ones. The time series <JA,x{t) corresponding to this 
scale is then the appropriate choice for further investiga¬ 
tions of the local dynamics. In our calculations we use 
circles with increasing radius for these areas. 

We describe now the procedure for the extraction of 
characteristic scales in four steps: 

i. Scatter plot of the time series {er„}: 

For the scatter plot we choose an embedding of 

Wn}- 

f(a n )=a n+k - (8) 

The lag k is determined by the condition that the 
auto-correlation of the time series is sufficiently 
small. This is fulfilled for k = 4. 

ii. Detecting candidates for unstable periodic orbits 
(UPOs): 


2 








x(n+k) 


The identification of unstable periodic orbits |[t|-|TT[ 
rests on the occurrence of a special sequence of 
points in the time series. If the system’s trajec¬ 
tory enters such a sequence, an UPO is visited by 
following a predictable pattern of values of the time 
series. An UPO is an intersection point of a stable 
and an unstable manifold Q|. In the vicinity of an 
UPO, we can approximate these manifolds locally 
linear. We define the following criteria for an UPO 
candidate (see Fig. [ilj): (1) the UPO itself is close to 
the line of identity: the perpendicular distance to 
the line of identity is smaller than the mean of the 
perpendicular distances of the five points and (2) 
a straight line approaches the line of identity (sta¬ 
ble manifold) and a straight line diverges from the 
line of identity (unstable manifold). We want to 
point out that these conditions are only capable to 
detect candidates for UPOs in time series, because 
they are necessary but not sufficient ones. We refer 
to jl^] for a method which is more appropriate to 
detect UPOs itself. 
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FIG. 2. Visual definition of an UPO candidate in a time 
series x(n)\ the dots are the data points, the squares denote 
the characteristic sequence of points (1,2,3: stable direction; 
3,4,5: unstable direction), point 3 is the UPO candidate which 
is close to the line of identity (dotted line). 


iii. Definition of the statistical significance: 

The statistical significance for the number of UPOs 
can be defined by comparing the original data with 
a large number of surrogate data files. Therefore, 
we count the number of UPOs in the original data, 
N, and in the surrogate data, N s , and compute the 
statistical significance [jTj: 


K = 


N - ( N s ) 
cr s 


(9) 


where a s is the standard deviation of N s and (N s ) 
the mean value of N s over all surrogate data files. 


For Gaussian distributions K > 3 is equivalent to a 
confidence level of 99% to reject the null hypothesis 
that the original data are linear in the sense that 
they do not contain a significant number of UPO 
candidates. Our surrogate data are generated by 
phase randomization and amplitude adjustment of 
the original data |Kj. This guarantees that the 
auto-correlation function as well as the distribution 
of the data are conserved. 

iv. Extraction of characteristic scales: 

The procedure to extract the characteristic scales 
from the data works as follows: One point in space 
is surrounded by circles with increasing radius: 
r = Ar, 2Ar, .... The step width is chosen to 
Ar = 5 km. For each circle a time series cr n is 
generated from the data as described in Sec. ||. 
In these time series we detect the UPO candidates 
and compute the significance K{r) in comparison 
with 100 surrogate data files. A characteristic scale 
around a point x with radius r max and signifi¬ 
cance K max is given, if (1) K{r ma x) > 3 and (2) 
K(r ma x) = Kmax = maxA'(r). This procedure 

r 

is applied to each point of a 20 x 20 lattice. To 
avoid finite size effects we exclude boundary points 
in space. Note also that the number of events in 
the boundary regions is too small to provide reli¬ 
able results. 

We have checked that the results do not depend sensi¬ 
tively on the parameters. 


IV. RESULTS OF DATA ANALYSIS 

Applying our algorithm to the earthquake data ex¬ 
hibits indeed in some regions such characteristic scales. 
Two examples are shown in Fig. pi, where we have plot¬ 
ted the significance K from Eq. (H) as a function of the 
spatial scale represented by an index for two fixed points 
in space. In Fig. ^(a) we observe values for K up to six 
and a clear maximum for the scale index 15. Here it is 
possible to assign a characteristic scale to the point x. In 
contrast to this we cannot find such a significant scale in 
Fig. Kb). The level of K = 3 is not reached here. Note 
that the significances in Fig. ||(a) and ||(b) are quite 
different, although the points x are very close in space. 
This result underlines that an averaging over the whole 
space or large parts of the space is indeed questionable. 

However, the rule for the extraction of the character¬ 
istic scale is rather simple. Not for all points a clear 
maximum of the significance can be observed. In some 
cases there are two or more maxima, or the shape of the 
curve K(r) is approximately constant. For the investi¬ 
gation of the local dynamics the function A should be 
studied in detail. But as a global approach, the rules for 
the selection of the scales seem to be reasonable. 
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FIG. 3. Significance K as a function of the spatial scale for 
two points in space: (a) x = (40.3,41.8), (b) x = (41.4,43.7). 
A unit of the scale index is 5 km. 

The distribution of the scales in space is given in Fig. [|. 
The scale sizes are related to the general seismic activ¬ 
ity (see Fig. [l]); the dynamics in very active regions is 
mostly determined on small scales, whereas in regions 
with less activity larger scales dominate. It is important 
to check, if this relationship is linear. In this case the 
distribution of the scales would simply reflect the seismic 
activity given in Fig. |l] and the results were completely 
trivial. This is, in fact, not the case, because the linear 
correlation coefficient between these quantities is almost 
zero. 

Our technique uses unstable periodic orbits to quan¬ 
tify nonlinear determinism. Next we check, whether a 
simpler discriminating statistic can be used for this aim. 
Therefore, we perform the same algorithm to obtain the 
time series (Eq. (||)), and compute then the skewness Jl(]] 
of these series instead of numbers of UPO candidates. 
The skewness is a rather simple quantity, which indi¬ 
cates nonlinear behavior and time reversal invariance in 
time series. Comparing original data and surrogate data 
by Eq. no significant deviation is observed. As a 
consequence, we need indeed a more powerful method to 


detect nonlinear features in the dynamics. 



42 44 46 48 

longitude 


FIG. 4. Distribution of characteristic scales in space 
(colour coded); a unit of the scale index is 5 km. 

V. CALCULATIONS WITH MODEL DATA 

To connect our results of data analysis with physical 
properties of the system, we want to apply our technique 
also to standard models of seismicity. In this way we can 
control the properties of the simulated data by changing 
some parameters of the model. 

An important point is the incompleteness of the earth¬ 
quake catalogue, which may have influence on the results. 
For instance the Gutenberg-Richter law is only fulfilled 
for 2 < to < 6 (see Fig. ||). In the region m > 6 the 
number of events is too small to provide a good statis¬ 
tics and for m < 2 a large number of micro-quakes is 
not measured for this catalogue. To analyse the depen¬ 
dence of the significance of UPOs on the completeness of 
the data, we generate a synthetic earthquake catalogue 
and assume that these surrogate earthquake data behave 
similar to the real data. 
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FIG. 5. The Gutenberg Richter law (Eq. (|l])): the points 
denote the number of earthquakes with magnitude greater or 
equal to m, the dashed line is a linear fit for 2 < m < 6. 

There is a standard model for earthquakes by Bak 
and Tang |l7j] exhibiting self-organized criticality. This 
is a cellular automaton analogue to well-known stick- 
slip-models |ll|, which describes the earth’s crust to 
be in a stationary critical state so that the distribution 
of earthquakes follows the Gutenberg-Richter relation. 
For a detailed description of the model we refer to 0- 
Here we only recall the rules for the redistribution of 
energy Z(i,j), if a cell ( i,j ) is in a critical state, i.e. 

^(bj) ^ ^crit '• 

Z(i,j)—>Z(i,j)~ 4 (10) 

Z(i±l,j±l)^Z(i±l,j±l) + l 

With a 50 x 50-grid and the critical value Z cr n = 3 we 
create two earthquake catalogues: (A) a “complete” cat¬ 
alogue including all micro-quakes and (B) a “truncated” 
one, where quakes with E < 10 (m < 1) are discarded. 
We choose the size of (2) equal to that of the real data. 
In Fig. H we compare the results of (A) and (B) for one 
point in space; like in Fig. ^ we compute the significance 
from Eq. (|)j) for UPOs as a function of the scale size. 
The shapes of the curves are very similar and the scales 
with hints for determinism are in almost all cases the 
same or at least very close to each other. Furthermore 
we see that this simple model yields scales with relatively 
high significances. While the dynamics for small scales 
is rather stochastic and provides only small significances, 
we observe the most deterministic behavior for interme¬ 
diate scales. For large scales the significance decreases 
again due to averaging effects in the dynamics, but does 
not tend to zero. 

However, in contrast to the real earthquake data the 
model is nearly homogeneous in space. For more realistic 
models one could add a component in Eq. (|To|) , which is 
space-dependent. 


FIG. 6. Significance K as a function of the spatial scale for 
the point x = (17,18) in the model of Bak and Tang; solid 
line: truncated model catalogue, dashed line: complete model 
catalogue. 

In the following we study another modification of the 
model, which yields a physical interpretation of the char¬ 
acteristic scales. Considering models with global cou¬ 
pling forces (e.g. slider block-models), the size of the 
scale should increase with growing coupling. To proof 
this for our model, we introduce a coupling strength C 
by changing the rules for the energy release. If Z(i,j) > 
Z C rit , the redistribution of the energy follows now 

Z(i,j)—> 0 (11) 

Z(i ± 1 ,j ± 1) — Z(i ± 1, j ± 1) + (1 - e)^Ml 

so that an increasing value of e corresponds to a decreas¬ 
ing coupling strength C := 1 — e. Note that Z is noninte¬ 
ger here in contrast to the model in |l7j. Moreover, a cell 
in a critical state decays by transferring its total energy 
to the neighbor cells. This modification is done in order 
to produce a large magnitude spectrum. In the model 
with integer energy units, the number of different magni¬ 
tudes decreases very fast with decreasing coupling. For 
different couplings we compute the total number of UPO 
candidates and compare it with 30 surrogate earthquake 
catalogues that are generated by randomizing the origi¬ 
nal catalogue such that the distribution is conserved. For 
one scale the mean number of UPO candidates is com¬ 
puted. In Fig. jrj(a) we show the statistical significance 
from Eq. (|j) for three different couplings C. We ob¬ 
serve a significant deviation between the model data and 
the mean of the surrogate data for intermediate scales. 
For small and large scales no significant determinism is 
present due to the statistics and averaging effects. Al¬ 
though the model is driven by stochastic forces, the dy¬ 
namics behaves deterministic on certain scales. More im¬ 
portant is the dependence of these scales on the coupling: 
for strong couplings we observe still positive significances 
on large scales; in contrast to this, the significance de¬ 
creases very fast in this scale range for small couplings. 
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On Fig. 0(b) the results for the calculation with the real 
data are given, which shows a qualitatively good agree¬ 
ment with the model with strong coupling in Fig. |(a). 




FIG. 7. (a) shows the significance for UPOs depending on 
the scale for three different couplings in the model described 
in the text (Eq. [n]) . For each scale the average number (over 
space) of UPOs is compared with 30 surrogate earthquake 
catalogues. The definition of the significance is analogue to 
Eq. (bf). One curve belongs to one catalogue, (b) results from 
the same calculation for the Armenia catalogue. A unit of the 
scale index is 10 km. 


Due to these results from pure model studies, we get 
indications that large spatial scales in real earthquake 
data correspond to strong coupling forces in the earth’s 
crust. 


eral approach for data analysis, because it is applicable 
to a large variety of systems. 

The main idea is to check whether it exists an interme¬ 
diate spatial scale between the noisy micro-scales and the 
large scales, where the dynamics is dominated by averag¬ 
ing effects. This intermediate scale is then characterized 
by a maximum of nontrivial determinism and it repre¬ 
sents the appropriate length scale, on which the main 
features of the underlying dynamics can be observed. 
To extract the characteristic spatial scale from the data, 
we look for unstable periodic orbits in time series cor¬ 
responding to different scales. The occurrence of such 
orbits is a measure for nonlinear and low-dimensional de¬ 
terminism. The scale with the highest significance (with 
respect to the condition K > 3) - in comparison with 
surrogate data - can be considered as the characteristic 
scale. 

Our calculations show that intermediate scales with 
deterministic dynamics in the sense as mentioned above 
exist in earthquake data. In some spatial regions we ob¬ 
serve a clear maximum for the significance as a function 
of the scale size. Moreover, in some regions the statistical 
significance reaches values up to seven. 

An interesting result arising from studies with model 
data is the interpretation of the spatial scales. Using a 
modification of the well-known model of Bak and Tang, 
we have shown that the scale is related to a coupling 
strength, i.e. high coupling strengths correspond to large 
scales. This seems to be an interesting tool to evaluate 
models, in particular such with locally varying coupling 
forces. 

However, one has to keep in mind that our ansatz is 
still very general and should be refined for special appli¬ 
cations. For instance some effects of seismicity are not 
yet included in our approach. Perhaps the technique can 
be improved by using more complicated areas than cir¬ 
cles. This would be well -adapted for the study of spatial 
inhomogeneities. In this context it is a special challenge 
to analyse the influence of such inhomogeneities in simple 
models. 

In the future we should focus on a detailed analysis 
of the time series. From seismology we know that every 
main shock is more or less accompanied by precursery 
phenomena |20] and aftershocks jl!|. It is an interesting 
question, if the occurrence of UPOs in the time series can 
be connected with these patterns. 

We believe, however, that our technique is promising 
for the analysis of a large variety of spatially extended 
natural systems. 


VI. SUMMARY AND OUTLOOK 

We have presented a new technique to characterize the 
dynamics of spatially extended natural systems. In par¬ 
ticular, we have analysed earthquake catalogues, but in 
principle the method can be regarded as a part of a gen- 
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